﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Text.RegularExpressions;
using System.Threading.Tasks;

namespace EarthGrid.core
{
    public class GridTool
    {
        #region 一些常量
        public const int xxx=256;
        #endregion
        /// <summary>
        /// 单位米
        /// </summary>
        public static int R_地球长半轴 = 6378137;

        /// <summary>
        /// 1弧度单位，1°对应的弧度为一个单位
        /// </summary>
        public static double R_1单位弧度 = Math.PI / 180;

        /// <summary>
        /// 前9位度的编码，然后6位分6位秒，后11位，共31位四进制
        /// </summary>
        public static int 最大网格编码长度 = 31;

        /// <summary>
        /// 108°58' 24.72"
        /// </summary>
        public static string 经纬度字符串正则 = @".*\d{1,3}[°,度]\s*\d{1,2}[',分]\s*\d+(\.\d+)?["", 秒].*";

        /// <summary>
        /// 度分秒格式的经纬度转换成十进制的经纬度
        /// </summary>
        /// <param name="str">必须包含度分秒的具体数值，格式与正则定义的匹配</param>
        /// <returns>十进制的经纬度值</returns>
        public static double 度分秒转十进制(string str)
        {
            if (Regex.IsMatch(str, 经纬度字符串正则))
            {
                string regStr = @"\d+(\.\d+)?";
                MatchCollection mc = Regex.Matches(str, regStr);
                Double d = 0D;
                int i = 0;
                int[] x = { 1, 60, 3600 };
                foreach (Match m in mc)
                {
                    Console.WriteLine(m.Value);
                    d += Double.Parse(m.Value) / x[i++];
                }
                return d;
            } 
            else
            {
                throw new Exception("请检查经纬度是否正确包含度分秒值");
            }
        }

        /// <summary>
        /// 根据有效经纬度判断实在哪个半球上
        /// </summary>
        /// <param name="lon">经度</param>
        /// <param name="lat">维度</param>
        /// <returns></returns>
        public static Hemisphere 判断在哪个半球(double lon, double lat)
        {
            if (lon > 0)
            {
                if (lat > 0)
                {
                    return Hemisphere.G0;
                }
                else
                {
                    return Hemisphere.G2;
                }
            }
            else
            {
                if (lat > 0)
                {
                    return Hemisphere.G1;
                }
                else
                {
                    return Hemisphere.G3;
                }
            }
        } 

        /// <summary>
        /// 从经纬度计算出所在精确网格
        /// </summary>
        /// <param name="lon"></param>
        /// <param name="lat"></param>
        /// <returns></returns>
        public static Grid 获取经纬度所在网格(LonLat lon, LonLat lat)
        {
            Int64 r = 经纬度编码交叉(lon, lat);
            
            String ra = 六十二位十进制转四进制(r);
            return null;
        }

        /// <summary>
        /// 东北半球上，定位精确度级网格编号
        /// </summary>
        /// <param name="lon"></param>
        /// <param name="lat"></param>
        /// <returns></returns>
        private static string G0GridFromDu(int lon, int lat)
        {
            return "";
        }

        /// <summary>
        /// 按规范要求，将经纬度进行莫顿交叉得到一个十进制数值
        /// </summary>
        /// <param name="lon"></param>
        /// <param name="lat"></param>
        /// <returns></returns>
        public static Int64 经纬度编码交叉(LonLat lon, LonLat lat)
        {

            // 先把度换成正十进制整数，秒后小数也转整数后转换成二进制取11位，那么小数后的最大值就是2047
            double lon_miao = Math.Round(lon.Miao, 4);
            int lon_xiaoshu = (int)((lon_miao - (int)lon_miao) * 1000);
            if (lon_xiaoshu > 2047)
            {
                lon_xiaoshu = 2047;
            }
            Int64 longitude = ((((((lon.Du) << 6) | lon.Fen) << 6) | (int)lon.Miao) << 11) | lon_xiaoshu;
            //Console.WriteLine("生成的31位二进制码：" + Convert.ToString(longitude, 2));
            double lat_miao = Math.Round(lat.Miao, 4);
            int lat_xiaoshu = (int)((lat_miao - (int)lat_miao) * 1000);
            Int64 latitude = ((((((lat.Du) << 6) | lat.Fen) << 6) | (int)lat.Miao) << 11) | (lat_xiaoshu > 2047 ? 2047 : lat_xiaoshu);
            //Console.WriteLine("生成的31位二进制码：" + Convert.ToString(latitude, 2));
            // 开始做错位交叉合并
            Int64 result = 0;
            for (int i = 0; i < 31; i++)
            {
                result |= ((longitude & (0b1 << i)) << i) | ((latitude & (0b1 << i)) << (i + 1));
               // Console.WriteLine(i+ " ： " + Convert.ToString(result, 2));
            }
            //Console.WriteLine("合并后的结果 ：  " + Convert.ToString(result, 2));
            return result;
        }

        /// <summary>
        /// 十进制转四进制，在这里不适应，编码长度会造成位数溢出
        /// </summary>
        /// <param name="a"></param>
        /// <returns></returns>
        public static Int64 十进制转四进制(int a)
        {
            //Console.WriteLine(Convert.ToString(a, 2));
            Int64 z = 0;
            int len = Convert.ToString(a, 2).Length;
            for (int i = 0; i <= len/2; i++)
            {
                int b = (int)a & (int)3;
                a >>= 2;
                z += b * (Int64)Math.Pow(10, i);
                /*Console.WriteLine("{0}, a : {1} , b : {2} , z : {3}", i, a, b, z);
                Console.WriteLine(Convert.ToString(a, 2));*/
            }
            return z;
        }

        /// <summary>
        /// 六十二位十进制转四进制
        /// </summary>
        /// <param name="a"></param>
        /// <returns></returns>
        public static String 六十二位十进制转四进制(Int64 a)
        {
            //Console.WriteLine(Convert.ToString(a, 2));
            int len = Convert.ToString(a, 2).Length / 2;
            String str = "";
            for (int i = 0; i < 31; i++)
            {
                int b = (int)a & (int)3;
                a >>= 2;
                str = b + str;
                /*Console.WriteLine("{0}, a : {1} , b : {2}, z : {3}", i, a, b, str);
                Console.WriteLine(Convert.ToString(a, 2));*/
            }
            return str;
        }

        /// <summary>
        /// 从已有网格分下级网格
        /// </summary>
        /// <param name="p"></param>
        /// <param name="立体">是否立体分割</param>
        /// <returns></returns>
        public static Grid[] createGridFromParent(Grid p, bool 立体)
        {
            if (p.Level >= 31)
            {
                return null;
            }
            
            // 定好下层四个网格数组的顺序，从最小坐标往最大坐标方向，先横向经度生成2个小格，再生成上面两小格
            Grid[] grids = p.InitNextLevel(false);

            (int LD_左下, int RD_右下, int LU_左上, int RU_右上) = getIndex(p.Hemisphere());

            grids[LD_左下].Code = p.Code + LD_左下;
            grids[RD_右下].Code = p.Code + RD_右下;
            grids[LU_左上].Code = p.Code + LU_左上;
            grids[RU_右上].Code = p.Code + RU_右上;

            // 判断是否到刚好到达分级或秒级网格，如果是，需要扩展60到64
            if (9 == p.Level || 15 == p.Level)
            {
                (Point p_LD, Point p_C, Point p_RU) = 扩展划分(p);
                grids[LD_左下].SetBox(p_LD, p_C);
                grids[RU_右上].SetBox(p_C, p_RU);
                grids[LU_左上].SetBox(p_LD.Lon, p_C.Lat, p_C.Lon, p_RU.Lat);
                grids[RD_右下].SetBox(p_C.Lon, p_LD.Lat, p_C.Lon, p_LD.Lat);
            }
            else
            {
                Point centerPoint = p.GetCenter();
                grids[LD_左下].SetBox(new Point(p.MinLon, p.MinLat), centerPoint);
                grids[RU_右上].SetBox(centerPoint, new Point(p.MaxLon, p.MaxLat));
                grids[LU_左上].SetBox(new Point(p.MinLon, centerPoint.Lat), new Point(centerPoint.Lon, p.MaxLat));
                grids[RD_右下].SetBox(centerPoint.Lon, p.MinLat, p.MaxLon, centerPoint.Lat);
            }
            if (立体)
            {
                // 高程方向进行分割
                return 高程方向分割(grids);
            }
            
            return grids;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="gs"></param>
        /// <returns></returns>
        private static Grid[] 高程方向分割(Grid[] gs)
        {
            Grid[] grids = new Grid[8];
            for (int i = 0; i < 4; i++)
            {
                grids[2 * i] = gs[i];
                grids[2 * i + 1] = gs[i].Copy();
                grids[2 * i].H_Code = gs[i].H_Code + "0";
                grids[2 * i + 1].H_Code = gs[i].H_Code + "1";
                if (gs[i].N_所在层 >= 0)
                {
                    grids[2 * i + 1].N_所在层 = gs[i].N_所在层 + 1;
                }
                else
                {
                    grids[2 * i + 1].N_所在层 = gs[i].N_所在层 - 1;
                }
            }
            return grids;
        }

        /// <summary>
        /// 根据网格所在半球，得出下级四个网格的数组下标
        /// </summary>
        /// <param name="hemisphere"></param>
        /// <returns></returns>
        private static (int i0, int i1, int i2, int i3) getIndex(Hemisphere hemisphere)
        {
            switch (hemisphere)
            {
                case Hemisphere.G0:
                    return (0, 1, 2, 3);
                case Hemisphere.G1:
                    return (1, 0, 3, 2);
                case Hemisphere.G2:
                    return (2, 3, 0, 1);
                case Hemisphere.G3:
                    return (3, 2, 1, 0);
                default:
                    return (0, 1, 2, 3);
            }
        }


        /// <summary>
        /// 网格在到分和秒级时，需要把60°扩展成64°，不同半球的扩展方向不同
        /// </summary>
        /// <param name="g"></param>
        /// <returns></returns>
        private static (Point pld, Point pcenter, Point pru) 扩展划分(Grid g)
        {
            Point pld = new Point();
            Point pcenter = new Point();
            Point pru = new Point();
            switch (g.Hemisphere())
            {
                case Hemisphere.G0:
                    pld.Lon = g.MinLon;
                    pld.Lat = g.MinLat;
                    if (9 == g.Level)
                    {
                        pcenter.Lon = new LonLat(LonLatEnum.lon, g.MinLon.Du, 32, 0);
                        pcenter.Lat = new LonLat(LonLatEnum.lat, g.MinLat.Du, 32, 0);
                        pru.Lon = new LonLat(LonLatEnum.lon, g.MinLon.Du, 64, 0);
                        pru.Lat = new LonLat(LonLatEnum.lat, g.MinLat.Du, 64, 0);
                    } 
                    else
                    {
                        pcenter.Lon = new LonLat(LonLatEnum.lon, g.MinLon.Du, g.MinLon.Fen, 32);
                        pcenter.Lat = new LonLat(LonLatEnum.lat, g.MinLat.Du, g.MinLat.Fen, 32);
                        pru.Lon = new LonLat(LonLatEnum.lon, g.MinLon.Du, g.MinLon.Fen, 64);
                        pru.Lat = new LonLat(LonLatEnum.lat, g.MinLat.Du, g.MinLat.Fen, 64);
                    }
                    break;
                case Hemisphere.G1:
                    if (9 == g.Level)
                    {
                        pld.Lon = new LonLat(LonLatEnum.lon, g.MaxLon.Du, -64, 0);
                        pld.Lat = g.MinLat;
                        pcenter.Lon = new LonLat(LonLatEnum.lon, g.MaxLon.Du, -32, 0);
                        pcenter.Lat = new LonLat(LonLatEnum.lat, g.MinLat.Du, 32, 0);
                        pru.Lon = g.MaxLon;
                        pru.Lat = new LonLat(LonLatEnum.lat, g.MinLat.Du, 64, 0);
                    }
                    else
                    {
                        pld.Lon = new LonLat(LonLatEnum.lon, g.MaxLon.Du, g.MaxLon.Fen, -64);
                        pld.Lat = g.MinLat;
                        pcenter.Lon = new LonLat(LonLatEnum.lon, g.MaxLon.Du, g.MaxLon.Fen, -32);
                        pcenter.Lat = new LonLat(LonLatEnum.lat, g.MinLat.Du, g.MinLat.Fen, 32);
                        pru.Lon = g.MaxLon;
                        pru.Lat = new LonLat(LonLatEnum.lat, g.MinLat.Du, g.MinLat.Fen, 64);
                    }
                    break;
                case Hemisphere.G2:
                    if (9 == g.Level)
                    {
                        pld.Lon = g.MinLon;
                        pld.Lat = new LonLat(LonLatEnum.lon, g.MaxLat.Du, -64, 0);
                        pcenter.Lon = new LonLat(LonLatEnum.lon, g.MinLon.Du, 32, 0);
                        pcenter.Lat = new LonLat(LonLatEnum.lat, g.MaxLat.Du, -32, 0);
                        pru.Lon = new LonLat(LonLatEnum.lon, g.MinLon.Du, 64, 0);
                        pru.Lat = g.MaxLat;
                    }
                    else
                    {
                        pld.Lon = g.MinLon;
                        pld.Lat = new LonLat(LonLatEnum.lon, g.MaxLat.Du, g.MaxLat.Fen, -64);
                        pcenter.Lon = new LonLat(LonLatEnum.lon, g.MinLon.Du, g.MinLon.Fen, 32);
                        pcenter.Lat = new LonLat(LonLatEnum.lat, g.MaxLat.Du, g.MaxLat.Fen, -32);
                        pru.Lon = new LonLat(LonLatEnum.lon, g.MinLon.Du, g.MinLon.Fen, 64);
                        pru.Lat = g.MaxLat;
                    }
                    break;
                case Hemisphere.G3:
                    if (9 == g.Level)
                    {
                        pld.Lon = new LonLat(LonLatEnum.lon, g.MaxLon.Du, -64, 0);
                        pld.Lat = new LonLat(LonLatEnum.lat, g.MaxLat.Du, -64, 0);
                        pcenter.Lon = new LonLat(LonLatEnum.lon, g.MaxLon.Du, -32, 0);
                        pcenter.Lat = new LonLat(LonLatEnum.lat, g.MaxLat.Du, -32, 0);
                    }
                    else
                    {
                        pld.Lon = new LonLat(LonLatEnum.lon, g.MaxLon.Du, g.MaxLon.Fen, -64);
                        pld.Lat = new LonLat(LonLatEnum.lat, g.MaxLat.Du, g.MaxLat.Fen, -64);
                        pcenter.Lon = new LonLat(LonLatEnum.lon, g.MaxLon.Du, g.MaxLon.Fen, -32);
                        pcenter.Lat = new LonLat(LonLatEnum.lat, g.MaxLat.Du, g.MaxLat.Fen, -32);
                    }
                    pru.Lon = g.MaxLon;
                    pru.Lat = g.MaxLat;
                    break;
            }
            return (pld, pcenter, pru);
        }

        #region 定义高程网格的公式

        /// <summary>
        /// 计算某层级网格距地心的距离
        /// </summary>
        /// <param name="i">网格所在层级</param>
        /// <param name="rad">网格的弧度</param>
        /// <returns></returns>
        public static double Ri_计算网格到地心距离(int i, double rad)
        {
            double r = 0;
            r = R_地球长半轴 * Math.Pow(1 + R_1单位弧度, i * rad/R_1单位弧度);
            return r;
        }

        /// <summary>
        /// 获取网格弧度，Π/180为一个弧度单位
        /// </summary>
        /// <param name="min">纬度起始值</param>
        /// <param name="max">纬度目标值</param>
        /// <returns></returns>
        public static double rad_弧度(double min, double max)
        {
            double len = Math.Abs(max - min);
            return len * Math.PI / 180;
        }

        /// <summary>
        /// 获取网格到地表距离
        /// </summary>
        /// <param name="g"></param>
        /// <returns></returns>
        public static double H_网格距地面高度(Grid g)
        {
            if (g.R_网格距地心距离 != 0)
            {
                return g.R_网格距地心距离 - R_地球长半轴;
            }
            double hudu = rad_弧度(g.MinLat.to十进制(), g.MaxLat.to十进制());
            double r = Ri_计算网格到地心距离(g.N_所在层, hudu);
            
            return  r - R_地球长半轴;
        }

        /// <summary>
        /// 计算网格的自身高
        /// </summary>
        /// <param name="g"></param>
        /// <returns></returns>
        public static double H_网格高度(Grid g)
        {
            double hudu = g.Rad_弧度;
            double h = Math.Pow((1 + hudu), g.N_所在层 * hudu / R_1单位弧度);
            h *= Math.Pow((1 + hudu), hudu / R_1单位弧度) - 1;
            return h;
        }
        #endregion

        #region 88纬度以上网格生成方法
        /// <summary>
        /// 生成第8级南北半球纬度88--90的两个网格
        /// </summary>
        public static Grid[] CreateLevel8GridLat88to90()
        {
            Grid gN = new Grid();
            Grid gS = new Grid();
            gN.Level = 8;
            gN.Code = "G00202200";
            gN.SetBox(new LonLat(LonLatEnum.lon, -180, 0, 0), new LonLat(LonLatEnum.lat, 88, 0, 0),
                    new LonLat(LonLatEnum.lon, 180, 0, 0), new LonLat(LonLatEnum.lat, 90, 0, 0));
            gS.Level = 8;
            gS.Code = "G20202200";
            gS.SetBox(new LonLat(LonLatEnum.lon, -180, 0, 0), new LonLat(LonLatEnum.lat, -90, 0, 0),
                    new LonLat(LonLatEnum.lon, 180, 0, 0), new LonLat(LonLatEnum.lat, -88, 0, 0));

            Grid[] gs = new Grid[2] { gN, gS };
            return gs;
        }

        /// <summary>
        /// 南北极中心圈的四分法，生成一个内圈，三个外层，这四个网格的分割方法不相同
        /// </summary>
        /// <param name="p0"></param>
        /// <returns>index 0 是内圈网格，2是跨正负经度的网格</returns>
        public static Grid[] 南北极四分法(Grid p0)
        {
            Grid[] gs = p0.InitNextLevel(false);
            // 极心位置 p00
            gs[0].Code = p0.Code + "2";
            // p01
            gs[1].Code = p0.Code + "0";
            // p02
            gs[2].Code = p0.Code + "1";
            // 判断是北半球还是南半球
            bool 北半球 = p0.MaxLat.Du > 0;

            // 外层经度三段，1（0° - 128°） 2（+128°-  -128°） 3（-128°- 0°）
            gs[0].MinLon = new LonLat(LonLatEnum.lon, -180, 0, 0);
            gs[0].MaxLon = new LonLat(LonLatEnum.lon, 180, 0, 0);

            gs[1].MinLon = new LonLat(LonLatEnum.lon, 0, 0, 0);
            gs[1].MaxLon = new LonLat(LonLatEnum.lon, 128, 0, 0);

            gs[2].MinLon = new LonLat(LonLatEnum.lon, 128, 0, 0);
            gs[2].MaxLon = new LonLat(LonLatEnum.lon, -128, 0, 0);

            gs[3].MinLon = new LonLat(LonLatEnum.lon, -128, 0, 0);
            gs[3].MaxLon = new LonLat(LonLatEnum.lon, 0, 0, 0);
            
            int center_fen = (p0.MinLat.Fen + p0.MaxLat.Fen) / 2;
            double center_miao = (p0.MinLat.Miao + p0.MaxLat.Miao) / 2;
            int center_du = (p0.MinLat.Du + p0.MaxLat.Du) / 2;

            if (北半球)
            {
                gs[3].Code = p0.Code.Replace("G0", "G1") + "0";
                if (p0.Level == 9)
                {
                    gs[0].MinLat = new LonLat(LonLatEnum.lat, p0.MinLat.Du, 32, 0);
                    gs[0].MaxLat = new LonLat(LonLatEnum.lat, p0.MinLat.Du, 64, 0);
                }
                else if (p0.Level == 15)
                {
                    gs[0].MinLat = new LonLat(LonLatEnum.lat, p0.MinLat.Du, p0.MinLat.Fen, -32);
                    gs[0].MaxLat = new LonLat(LonLatEnum.lat, p0.MinLat.Du, p0.MinLat.Fen, -64);
                }
                else
                {
                    gs[0].MinLat = new LonLat(LonLatEnum.lat, center_du, center_fen, center_miao);
                    gs[0].MaxLat = p0.MaxLat;
                }
                gs[1].MinLat = gs[2].MinLat = gs[3].MinLat = p0.MinLat;
                gs[1].MaxLat = gs[2].MaxLat = gs[3].MaxLat = gs[0].MinLat;
            }
            else
            {
                gs[3].Code = p0.Code.Replace("G2", "G3") + "0";
                if (p0.Level == 9)
                {
                    gs[0].MinLat = new LonLat(LonLatEnum.lat, p0.MaxLat.Du, -64, 0);
                    gs[0].MaxLat = new LonLat(LonLatEnum.lat, p0.MaxLon.Du, -32, 0);
                }
                else if (p0.Level == 15)
                {
                    gs[0].MinLat = new LonLat(LonLatEnum.lat, p0.MaxLat.Du, p0.MaxLat.Fen, -64);
                    gs[0].MaxLat = new LonLat(LonLatEnum.lat, p0.MaxLon.Du, p0.MaxLat.Fen, -32);
                }
                else
                {
                    gs[0].MinLat = p0.MinLat;
                    gs[0].MaxLat = new LonLat(LonLatEnum.lat, center_du, center_fen, center_miao);
                }
                gs[1].MinLat = gs[2].MinLat = gs[3].MinLat = gs[0].MaxLat;
                gs[1].MaxLat = gs[2].MaxLat = gs[3].MaxLat = p0.MaxLat;
            }
            return gs;
        }

        /// <summary>
        /// 南北极圈四分后的P2的四分法
        /// </summary>
        /// <param name="p2"></param>
        /// <returns></returns>
        public static Grid[] 南北极P2跨正负的四分法(Grid p2)
        {
            Grid[] gs = p2.InitNextLevel(false);
            bool 北半球 = p2.MaxLat.Du > 0;
            // 按经度分实际应该是 (128°  -  180°)  (-180°  -  -128°)
            // 为方便等分计算，东经两个应扩展到(128°  -  192°)  (-192° -   -128°)
            gs[0].MinLon = new LonLat(LonLatEnum.lon, 128, 0, 0);
            gs[0].MaxLon = new LonLat(LonLatEnum.lon, 192, 0, 0);
            gs[0].Code = p2.Code + "0";
            gs[1].MinLon = gs[0].MinLon;
            gs[1].MaxLon = gs[0].MaxLon;
            gs[1].Code = p2.Code + "2";

            // 西经两个
            string cp = 北半球 ? p2.Code.Replace("G0", "G1") : p2.Code.Replace("G1", "G3");
            gs[2].MinLon = new LonLat(LonLatEnum.lon, -192, 0, 0);
            gs[2].MaxLon = new LonLat(LonLatEnum.lon, -128, 0, 0);
            gs[2].Code = cp + "0";
            gs[3].MinLon = gs[2].MinLon;
            gs[3].MaxLon = gs[2].MaxLon;
            gs[3].Code = cp + "2";

            // 0,2的纬度低 1,3的纬度高
            if (北半球)
            {
                
                if (p2.Level == 9)
                {
                    gs[0].MaxLat = gs[2].MaxLat = new LonLat(LonLatEnum.lat, p2.MinLat.Du, 32, 0);
                    gs[1].MaxLat = gs[3].MaxLat = new LonLat(LonLatEnum.lat, p2.MinLat.Du, 64, 0);
                }
                else if (p2.Level == 15)
                {
                    gs[0].MaxLat = gs[2].MaxLat = new LonLat(LonLatEnum.lat, p2.MinLat.Du, p2.MinLat.Fen, 32);
                    gs[1].MaxLat = gs[3].MaxLat = new LonLat(LonLatEnum.lat, p2.MinLat.Du, p2.MinLat.Fen, 64);
                }
                else
                {
                    gs[0].MaxLat = gs[2].MaxLat = new LonLat(LonLatEnum.lat, p2.MinLat.Du,
                        (p2.MinLat.Fen + p2.MaxLat.Fen) / 2, (p2.MinLat.Miao + p2.MaxLat.Miao) / 2);
                    gs[1].MaxLat = gs[3].MaxLat = p2.MaxLat;
                }
                gs[0].MinLat = gs[2].MinLat = p2.MinLat;
                gs[1].MinLat = gs[3].MinLat = gs[0].MaxLat;
            }
            else
            {
                if (p2.Level == 9)
                {
                    gs[0].MinLat = gs[2].MinLat = new LonLat(LonLatEnum.lat, p2.MaxLat.Du, -32, 0);
                    gs[1].MinLat = gs[3].MinLat = new LonLat(LonLatEnum.lat, p2.MaxLat.Du, -64, 0);
                }
                else if (p2.Level == 15)
                {
                    gs[0].MinLat = gs[2].MinLat = new LonLat(LonLatEnum.lat, p2.MaxLat.Du, p2.MaxLat.Fen, -32);
                    gs[1].MinLat = gs[3].MinLat = new LonLat(LonLatEnum.lat, p2.MaxLat.Du, p2.MaxLat.Fen, -64);
                }
                else
                {
                    gs[0].MinLat = gs[2].MinLat = new LonLat(LonLatEnum.lat, p2.MaxLat.Du,
                        (p2.MaxLat.Fen + p2.MinLat.Fen) / 2, (p2.MaxLat.Miao + p2.MinLat.Miao) / 2);
                    gs[1].MinLat = gs[3].MinLat = p2.MinLat;
                }
                gs[0].MaxLat = gs[2].MaxLat = p2.MaxLat;
                gs[1].MaxLat = gs[3].MaxLat = gs[0].MinLat;
            }
            return gs;
        }

        public static Grid[] 南北极P1P3的四分法(Grid p)
        {
            Grid[] gs = p.InitNextLevel(false);
            (LonLat l_lat, LonLat c_lat, LonLat h_lat) = 纬度分割(p);
            LonLat c_lon = new LonLat(LonLatEnum.lon, (p.MinLon.Du + p.MaxLon.Du) / 2, 0, 0); // TODO 经度的分法需要重新考虑
            
            gs[0].Code = p.Code + "0";
            gs[1].Code = p.Code + "1";
            gs[2].Code = p.Code + "2";
            gs[3].Code = p.Code + "3";

            switch (p.Hemisphere())
            {
                case Hemisphere.G0:
                    gs[0].SetBox(p.MinLon, l_lat, c_lon, c_lat);
                    gs[1].SetBox(c_lon, l_lat, p.MaxLon, c_lat);
                    gs[2].SetBox(p.MinLon, c_lat, c_lon, h_lat);
                    gs[3].SetBox(c_lon, c_lat, p.MaxLon, h_lat);
                    break;
                case Hemisphere.G1:
                    gs[0].SetBox(c_lon, l_lat, p.MaxLon, c_lat);
                    gs[1].SetBox(p.MinLon, l_lat, c_lon, c_lat);
                    gs[2].SetBox(c_lon, c_lat, p.MaxLon, h_lat);
                    gs[3].SetBox(p.MinLon, c_lat, c_lon, h_lat);
                    break;
                case Hemisphere.G2:
                    gs[0].SetBox(p.MinLon, c_lat, c_lon, l_lat);
                    gs[1].SetBox(c_lon, c_lat, p.MaxLon, l_lat);
                    gs[2].SetBox(p.MinLon, h_lat, c_lon, c_lat);
                    gs[3].SetBox(c_lon, h_lat, p.MaxLon, c_lat);
                    break;
                case Hemisphere.G3:
                    gs[0].SetBox(c_lon, c_lat, p.MaxLon, l_lat);
                    gs[1].SetBox(p.MinLon, c_lat, c_lon, l_lat);
                    gs[2].SetBox(c_lon, h_lat, p.MaxLon, c_lat);
                    gs[3].SetBox(p.MinLon, h_lat, c_lon, c_lat);
                    break;
            }
            return gs;
        }

        /// <summary>
        /// 南北极第10级以后的外层网格分法，第9级的三个外层网格分割时，调用特定的分法后得到第10级网格，
        /// 然后就可以使用这个方法来分那些10级网格
        /// 第10级的经度跨度是64°，到16级时经度跨度是1°，17级时跨度是32′，22级是1′，23级是32″，秒以后没有扩展
        /// </summary>
        /// <param name="g"></param>
        /// <returns></returns>
        public static Grid[] 南北极10级以下外层分法(Grid g)
        {
            Grid[] gs = g.InitNextLevel(false);
            (LonLat l_lat, LonLat c_lat, LonLat h_lat) = 纬度分割(g);
            (LonLat min_lon, LonLat center_lon, LonLat max_lon) = 经度分割(g);
            gs[0].Code = g.Code + "0";
            gs[1].Code = g.Code + "1";
            gs[2].Code = g.Code + "2";
            gs[3].Code = g.Code + "3";
            switch (g.Hemisphere())
            {
                case Hemisphere.G0:
                    gs[0].SetBox(min_lon, l_lat, center_lon, c_lat);
                    gs[1].SetBox(center_lon, l_lat, max_lon, c_lat);
                    gs[2].SetBox(min_lon, c_lat, center_lon, h_lat);
                    gs[3].SetBox(center_lon, c_lat, max_lon, h_lat);
                    break;
                case Hemisphere.G1:
                    gs[0].SetBox(center_lon, l_lat, max_lon, c_lat);
                    gs[1].SetBox(min_lon, l_lat, center_lon, c_lat);
                    gs[2].SetBox(center_lon, c_lat, max_lon, h_lat);
                    gs[3].SetBox(min_lon, c_lat, center_lon, h_lat);
                    break;
                case Hemisphere.G2:
                    gs[0].SetBox(min_lon, c_lat, center_lon, l_lat);
                    gs[1].SetBox(center_lon, c_lat, max_lon, l_lat);
                    gs[2].SetBox(min_lon, h_lat, center_lon, c_lat);
                    gs[3].SetBox(center_lon, h_lat, max_lon, c_lat);
                    break;
                case Hemisphere.G3:
                    gs[0].SetBox(center_lon, c_lat, max_lon, l_lat);
                    gs[1].SetBox(min_lon, c_lat, center_lon, l_lat);
                    gs[2].SetBox(center_lon, h_lat, max_lon, c_lat);
                    gs[3].SetBox(min_lon, h_lat, center_lon, c_lat);
                    break;
            } 
            return gs;
        }

        /// <summary>
        /// 南北极纬度的分割，南北纬89~90度之间
        /// </summary>
        /// <param name="g"></param>
        /// <returns>返回中低高三个纬度值</returns>
        private static (LonLat l_lat, LonLat c_lat, LonLat h_lat) 纬度分割(Grid g)
        {
            // 低纬度
            LonLat l_lat;
            // 中纬度
            LonLat c_lat;
            // 高纬度
            LonLat h_lat;
            if (g.MaxLat.Du > 0)
            {
                l_lat = g.MinLat;
                if (g.Level == 9) 
                {
                    c_lat = new LonLat(LonLatEnum.lat, l_lat.Du, 32, 0);
                    h_lat = new LonLat(LonLatEnum.lat, l_lat.Du, 64, 0);
                }
                else if (g.Level == 15)
                {
                    c_lat = new LonLat(LonLatEnum.lat, l_lat.Du, l_lat.Fen, 32);
                    h_lat = new LonLat(LonLatEnum.lat, l_lat.Du, l_lat.Fen, 64);
                }
                else
                {
                    c_lat = new LonLat(LonLatEnum.lat, l_lat.Du, (l_lat.Fen + g.MaxLat.Fen) / 2, (l_lat.Miao + g.MaxLat.Miao) / 2);
                    h_lat = g.MaxLat;
                }
            }
            else
            {
                l_lat = g.MaxLat;
                if (g.Level == 9)
                {
                    c_lat = new LonLat(LonLatEnum.lat, l_lat.Du, -32, 0);
                    h_lat = new LonLat(LonLatEnum.lat, l_lat.Du, -64, 0);
                }
                else if (g.Level == 15)
                {
                    c_lat = new LonLat(LonLatEnum.lat, l_lat.Du, l_lat.Fen, -32);
                    h_lat = new LonLat(LonLatEnum.lat, l_lat.Du, l_lat.Fen, -64);
                }
                else
                {
                    c_lat = new LonLat(LonLatEnum.lat, l_lat.Du, (g.MinLat.Fen + g.MaxLat.Fen) / 2, (g.MinLat.Miao + g.MaxLat.Miao) / 2);
                    h_lat = g.MinLat;
                }
            }
            return (l_lat, c_lat, h_lat);
        }

        /// <summary>
        /// 南北极纬度在88~90度时，经度合并上的分割
        /// </summary>
        /// <param name="g"></param>
        /// <returns></returns>
        private static (LonLat min_lon, LonLat center_lon, LonLat max_lon) 经度分割(Grid g)
        {
            LonLat min_lon;
            LonLat center_lon;
            LonLat max_lon;

            if (g.Level == 16)
            {
                min_lon = g.MaxLon.Du > 0 ? g.MinLon : new LonLat(LonLatEnum.lon, g.MaxLon.Du, -64, 0);
                center_lon = g.MaxLon.Du > 0 ? new LonLat(LonLatEnum.lon, g.MinLon.Du, 32, 0) : new LonLat(LonLatEnum.lon, g.MaxLon.Du, -32, 0);
                max_lon = g.MaxLon.Du > 0 ? new LonLat(LonLatEnum.lon, g.MinLon.Du, 64, 0) : new LonLat(LonLatEnum.lon, g.MaxLon.Du, -64, 0);
            } 
            else if (g.Level == 22)
            {
                min_lon = g.MaxLon.Du > 0 ? g.MinLon : new LonLat(LonLatEnum.lon, g.MaxLon.Du, g.MaxLon.Fen, -64);
                center_lon = g.MaxLon.Du > 0 ? new LonLat(LonLatEnum.lon, g.MinLon.Du, g.MinLon.Fen, 32) : new LonLat(LonLatEnum.lon, g.MaxLon.Du, g.MaxLon.Fen, -32);
                max_lon = g.MaxLon.Du > 0 ? new LonLat(LonLatEnum.lon, g.MinLon.Du, g.MinLon.Fen, 64) : new LonLat(LonLatEnum.lon, g.MaxLon.Du, g.MaxLon.Fen, -64);
            } 
            else
            {
                min_lon = g.MinLon;
                max_lon = g.MaxLon;
                center_lon = new LonLat(LonLatEnum.lon, (max_lon.Du + min_lon.Du) / 2, (min_lon.Fen + max_lon.Fen) / 2, (min_lon.Miao + max_lon.Miao) / 2);
            }
            
            return (min_lon, center_lon, max_lon);
        }

        #endregion
    }
}
